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We theoretically investigate the interaction between a single molecule and a metallic nanoparticle. 
We develop a general quantum mechanical description for the calculation of the enhancement of 
radiative and non-radiative decay channels for a molecule situated in the nearfield regime of the 
metallic nanoparticle. Using a boundary element method approach, we compute the scattering 
rates for several nanoparticle shapes. We also introduce an eigenmode expansion and quantization 
scheme for the surface plasmons, which allows us to analyze the scattering processes in simple 
physical terms. An intuitive explanation is given for the large quantum yield of quasi one- and two- 
dimensional nanostructures. Finally, we briefly discuss resonant Forster energy transfer in presence 
of metallic nanoparticles. 
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I. INTRODUCTION 

Metal nanoparticles can sustain local surface plasmon 
excitations, particle plasmons, which are hybrid modes of 
a light field coupled to a coherent electron charge oscil- 
lation [l] [2] . The properties of these excitations depend 
strongly on particle geometry and interparticle coupling, 
and give rise to a variety of effects, such as frequency- 
dependent absorption and scattering or near field en- 
hancement. Particle plasmons enable the concentration 
of light fields to nanoscale volumes and play a key role in 
surface enhanced spectroscopy [3 . This can be exploited 
for a variety of applications, such as surface-enhanced 
Raman scattering [4, 5 or biochemical sensorics [6 . 

Improved nanofabrication methods nowadays allow ad- 
vanced control of the nanoparticle shape and the arrange- 
ment of nanoparticle ensembles [7 , and open the possi- 
bility to flexibly tailor specific molecule-nanoparticle cou- 
plings. Anger et al. [8^ recently investigated the fluores- 
cence of a single molecule and a single laser-irradiated 
spherical gold nanoparticle, and demonstrated by vary- 
ing the molecule-nanoparticle distance the continuous 
transition from fluorescence enhancement to fluorescence 
quenching [9^. This finding was supported by Kiihn and 
coworkers [10 using a similar setup, and stimulated a 
large number of related experiments [lTJ[T2]. 

The theoretical description of the modified molecule 
dynamics in presence of metallic nanoparticles (MNPs) 
has been addressed by several authors in the past (see, 
e.g. [T, HH [Ei and references therein) . More recently, 
the question of how to treat quantum electrodynamics in 
the presence of absorbing media has seen anewed inter- 
est [ini [T71 [iSl |19] , motivated by related work on thermal 
nearfield radiation [20l[2T] and thermally induced scatter- 
ing processes of atoms located in the vicinity of metallic 
or superconducting bodies |22l [23l [24l [25]. This work 
has identified the Green tensor as the central element to 
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establish, in linear response, a convenient link between 
quantum electrodynamics and classical Maxwell theory. 

In this work we start by reviewing the description 
scheme for a molecular dipole placed in the vicinity of 
a metallic nanoparticle. We show how to use the dyadic 
Green function of Maxwell's theory to express the elec- 
tric field fluctuations in terms of the current fluctuations 
in the metal, and how to derive the molecule scatter- 
ing rate within linear response theory. In the quasistatic 
regime, where the light wavelength is much larger than 
the nanoparticle size and the molecule-MNP distance, 
the decay rate can be separated into a radiative and non- 
radiative contribution, which both show a simple scal- 
ing behavior. We next use a boundary element method 
approach [26 to compute these decay rates for realistic 
nanoparticle geometries. The results are shown to be 
in nice agreement with Mie theory. We additionally in- 
troduce for a Drude description of the metal electrons 
a quantization scheme for the surface plasmons [271 HI] • 
Although this approach could be used to study the non- 
linear optical response of metallic nanoparticles, in this 
work we use this framework to provide an intuitive pic- 
ture for the radiative and non-radiative molecular de- 
cay channels only, and give a simple explanation for the 
higher quantum yield achievable in quasi one- and two- 
dimensional nanoparticles. Finally, we briefly discuss res- 
onant Forster energy transfer between donor and accep- 
tor molecules placed in the nanoparticle nearfield. 



II. THEORY 

As for the description of the molecule, in this work 
we follow Refs. [Mj [15] who considered a generic two- 
level system. This approach is also best suited for other 
quantum emitters, such as colloidal quantum dots. We 
describe the molecule in terms of a generic two level 
scheme, with the groundstate and 1 the excited molecu- 
lar state with energy Ei. Other, more refined description 
schemes, considering, e.g., the different vibronic states of 
the molecule [29] [30] or details of the molecular orbitals 
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[31], could be treated in a similar fashion. 

We first consider the optical decay from the excited 
state 1 to the groundstate. Associated to this transition 
is a dipole operator 

d = (a_ + a+) = d+ + d" , (1) 

where /j. is the molecular dipole moment, assumed to 
be real, (j+ = |1)(0| is the molecular excitation operator 
which brings the molecule from the ground to the excited 
state, and &- = |0)(1| the corresponding de-excitation 
operator. In an interaction representation according to 
the free Hamiltonian = the operators a_ 

and (j+ oscillate with positive and negative frequencies 
e"*'^^ and e*'^^, respectively, as indicated by the dipole 
de-excitation and excitation operators in Eq. ([T]). It 
is clear that the latter operators can be easily generalized 
in case of more complicated level schemes. 

In a similar fashion, we split the electric field operator 
E{r) = E^ir) + E~{r) into two contributions evolv- 
ing with positive and negative frequencies [32^^ "33^ . The 
molecule-light coupling is described within the dipole ap- 
proximation 

Hop = -d-E^- (d-^ -E- ^d- ■ E^^ , (2) 

where E{rm) has to be taken at the position Vm of the 
molecule. We have neglected contributions where both d 
and E oscillate with either positive or negative frequen- 
cies, which corresponds to the common rotating-wave 
approximation [32l [33] . In quantum optics one usually 
expands the field operators E^ in the eigenmodes of 
Maxwell's equation to obtain photon annihilation and 
creation operators. The first expression in parentheses 
then describes a molecular de-excitation and the creation 
of a photon, and the second expression the reversed pro- 
cess. 

Things change considerably if the field operators E 
evolve in a complex dielectric environment where absorp- 
tion takes place, as is the case for MNPs. Owing to ab- 
sorption the eigenmodes acquire a complex energy, the 
imaginary part associated with the damping of the eigen- 
modes, which spoils the usual quantization procedure for 
the electromagnetic fields. As we will show next, in linear 
response one can make use of the fluctuation- dissipation 
theorem ^ to avoid such an eigenmode expansion of E. 



A. Molecule-MNP coupling 

1. Green function approach 

Our starting point is provided by the Fermi-golden-rule 
expression for the molecule decay rate 

7 = 27r/ii (El{rm,(^)EJ{rm,(^)) l-ij (3) 



where we have made use of the Einstein sum convention.^ 
Some details of its derivation are provided in appendix [A] 
The expression in brackets describes the field fluctuations 
in thermal equilibrium, to be determined at the molecu- 
lar transition frequency ix) and position r^. These fluc- 
tuations are induced by current noise in the metal. Let 
jf(r) = j^{r) -\-j~{r) denote the current operator in the 
metal, which we have again split into positive and nega- 
tive frequency components. If j were a classical current, 
evolving with frequency e~*^^ the electric field could be 
determined from the wave equation [35 

V X V X E{r, uj) - k^e{r, uj)E{r, uj) = j(r, uj) , (4) 

where E{r^uo) and j{r^u;) are the Fourier transform of 
the electric field and current, respectively, c is the vac- 
uum speed of light, k = uo/c the wavenumber in vacuum, 
and e{r^uo) the space and frequency dependent dielec- 
tric function. We set the magnetic permeability ji = 1 
throughout. A convenient way to compute the electric 
field for a given current source is by means of the dyadic 
Green tensor G{r^r'^uo) [2 , to be determined from the 
wave equation 

V X V X G(r, r'uo) - k^e{r, cj)G(r, r' , uo) = 47rJ(r -r')t, 

(5) 

together with appropriate boundary conditions. The cal- 
culation of G(r,r',a;) is a common problem in classical 
electrodynamics, and will be described in more detail fur- 
ther below. For the moment it suffices to assume that the 
Green function is at our hand. By comparison of Eqs. (|4| 
and ([5| we readily observe 

E{r, ^) = 5 / G{r, r'Lo)j{r', lo) dr' . (6) 

Thus, for a linear material response, described by the 
dielectric function e(r,cc;), the Green function provides 
the link between the current source and the electric field. 
Eq. (|6| can be also used for the corresponding electric- 
field and current operators [16 , where E^{r^u) is in- 
duced by j^{r\ij). Note that in the linear response 
regime positive and negative frequency components do 
not mix. We can now use Eq. ([6| to relate the field fluc- 
tuations through 

X (J + (r,c^)j^(r^c^))G*,(r^,r^c^)^r^r' (7) 

to the current fluctuations in the dielectric. Here, the 
Green functions describe how the field propagates from 
the current source to the position of the molecule. In 
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appendix [B] we show that for local and isotropic di- 
electrics the current correlation function can be expressed 
as SkiS{r — r')uj'^e"{r,u)/(2TT), with e"{ijj) being the 
imaginary part of the dielectric function. Together with 
the integral equation [16^ 



Gik(r, s, uj) k e\s, uj) G]^{r' , s, cj) ds 

= A^^m[G^j{r,r',u)] , (8) 

which is obtained by multiplying Eq. ([5| with G*^ and 
subtracting the equation obtained by complex conjuga- 
tion, we then arrive at the expression for the molecular 
decay rate 



(9) 



Equation ([9| is our final result. It shows that the decay 
rate of a molecule in proximity to a MNP is fully deter- 
mined by the dyadic Green function of classical Maxwell 
theory. This is a huge simplification because all the dy- 
namics of current fiuctuations in the metal are embodied 
in the dielectric function, which can be obtained from 
either experiment or first-principles calculation. 

Eq. ([9| can be brought into an even more transparent 
form by using the relation j = —iufi between the current 
and the molecular dipole. Inserting this expression into 
Eq. (|6|, one immediately obtains for Eq. (|9| the form 



7 
2 



(10) 



where E{rm) is the electric field induced by the dipole ii 
itself. Equation (10) describes a self-interaction, where 



the dipole polarizes the MNP, and the total electric field 
of the dipole and of the polarized MNP acts, in turn, 
back on the dipole. It is the imaginary part of this self 
interaction which accounts for the decay of the excited 
molecular state. 



molecular dipole, as computed in the quasistatic limit. 
We have to additionally account for radiative losses. This 
is done by using the standard expression 
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+ MmnpI 



(12) 



for the composite molecule-MNP dipole radiator. Here 
7o = |n5/i^/c^ is the free-space decay rate of the molecule, 
with ni) being the refractive index of the embedding 
medium. The second term on the right-hand side of 



Eq. (12) is the enhancement of the radiative molecular 



decay in presence of the MNP, where /l^mnp is the MNP 
dipole induced by the molecule. Below we will show how 
to compute £^qs(r^) and /l^mnp within a boundary ele- 
ment method approach. 



3. Scaling 

In the general case there are two length scales which 
define the problem, namely the size of the nanoparticle 
and the wavelength of the light. In the quasistatic regime, 
where the light wavelength is assumed to be much larger 
than the MNP size, it is a single length scale entering the 
problem. In this regime one can find a simple scaling law 
for 7nr and 7r on the size of the MNP. Suppose that we 
have computed the two decay rates for a given nanopar- 
ticle geometry. We can then scale the whole geometry by 
a scaling factor A, according to r ^ Ar, while keeping all 
other quantities, such as transition frequency uo or molec- 
ular dipole moment /l^, fixed. One can show that such a 
scaling transforms the scattering rates according to 



7nr(A) = 7nr/A^ , 7r(A) = 7r 



(13) 



2. Quasistatic approximation 

In many cases of interest the size of the MNP is much 
smaller than the light wavelength. It is then possible 
to employ the so-called quasistatic approximation. For 
a given source, the electric field is computed from the 
scalar potential ^{r) according to E{r) = — V<I>(r). This 
approach differs from the truly static case in that one 
uses the frequency-dependent dielectric function e(r,cj), 
at the molecular transition frequency cj, rather than the 
static limit. The imaginary part of e(r,a;) is associated 
with Ohmic losses of electromagnetic fields inside the 
metal. The self-interaction expression in the quasistatic 
limit thus accounts only for nonradiative losses 



We will find that Eq. (13) is an extremely useful ex- 



7n 



(11) 



as indicated by the subscript on 7. £^qs(r^) is the elec- 
tric field at the position of the molecule, induced by the 



pression for estimating the relative importance of non- 
radiative and radiative decay rates for realistic MNPs. 



B. Mie theory 

The electromagnetic response of spherical metallic par- 
ticles can be solved exactly within Mie theory. This the- 
ory provides closed expressions in both the retarded and 
quasistatic case. Indeed, it is a very rare situation that 
one can find an exact solution for an experimentally inter- 
esting, nontrivial problem. In this work we will be using 
the Mie theory primarily as a test case for our numerical 
scheme described below. As there exists a vast amount 
of literature on Mie theory [36 , in appendix [C] we only 
briefiy sketch the quasistatic approach and provide the 
explicit expressions used in our calculations. 
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C. Boundary element method (BEM) 

In the quasistatic hmit the wave equation Q for the 
electric field transforms to the Poisson equation 



V^^(r) = -47rp(r) 



(14) 



for the scalar potential. Here p{r) is the charge distri- 
bution of the molecular dipole. The problem can be sig- 
nificantly simplified for a MNP described by a homoge- 
neous dielectric function em{^) embedded in a medium 
with dielectric constant e^. As long as there is no danger 
of confusion, we will suppress the frequency dependence 
of the metal dielectric function and assume that is 
evaluated at the molecular transition frequency We can 
now use the static Green function G{r^r') = l/(|r — r'l), 
which is the solution of the Poisson equation 



where the different signs depend on whether the surface 
derivative is taken inside or outside the nanoparticle. 

At this point it is convenient to change from boundary 
integrals to boundary elements, suitable for a compu- 
tational implementation. Within the boundary element 
method (BEM) approach, the surface of the nanoparticle 
is approximated by small surface elements of triangular 
shape, as shown in Fig. 1 for the example of a spheri- 
cal particle. In its most simple form, which we will be 
using in this work, the surface charges are assumed to 
be located at the center of each triangle s^, and the ma- 
trices Gij and Fij connecting different surface elements 
are obtained from G{si^Sj) and F{si^Sj).^ The integral 
equation (17) then reduces to two matrix equations 



$:„ = (F + 27rl)(T + $; 



ext , *b = iF- 27rl)a ^ , 



y^G{r,r') = -4n5{r-r') 



(15) 



for a point source located at position r^ The solution of 
Eq. ( 14 ) can then be expressed as [26l |37] 



(18) 

where F is the matrix with elements F{si^Sj) and the 
subscript of indicates the side on which the surface 
derivative is taken. We can now use the boundary con- 
dition em^'rn ~ ^b^b = to computc, for a given external 
charge distribution, the surface charges according to 



where dCt denotes the surface of the metallic nanopar- 
ticle, cr is a surface charge that has to be chosen such 
that the appropriate boundary conditions are fulfilled, 
and $ext is the scalar potential for the external charge 
distribution of the molecular dipole. Indeed, Eq. (16) ful- 



Equation ( |19[ ) is the main result for our BEM ap- 
proach. It shows that the surface charge can be com- 
puted through numerical inversion of a matrix having the 
size of the number of triangular surface elements, which 
is typically of the order of thousands. Once we have 
fills the Poisson equation ^ everywhere excep^at the computed a for a given inhomogeneity we can com- 
pute the potential and electric fields at any given space 
point from Eq. ( 16 ), and the induced dipole moment from 
Mmnp = ^i^i' computational approach we de- 

scribe the molecular dipole by two oppositely charged 
particles separated by a small distance, and compute the 
electric field from Eq. ( 16 ) by means of a finite difference 
scheme. 



boundaries. The boundary conditions at dQ are the con- 
tinuity of the tangential electric field and h-{Di) — Dm) = 
0, where n is the unit vector normal to the surface and 
pointing away from the MNP. is the dielectric dis- 
placement at the outer surface, and Dm the displacement 
at the inner surface. The continuity of the electric field 
is guaranteed for a potential which is continuous at the 
boundary, as is the case in Eq. (16). To implement the 



second boundary condition, we have to take the surface 
derivative of $ at either side of the boundary. The limit 
r ^ s in Eq. (16) has to be treated with some care.^ 
With the abbreviations F(s, s') = (n • V)G(s, s') and 
= (n • V)^, we then get 

lim<I>'(r)= / F{s,s')a{s')ds' ±27Ta{s) ^^'^^^{s) , 

JdQ 

(17) 



D. Eigenmode expansion 

1. Boundary element method approach 

Things can be formulated differently when the dielec- 
tric function is approximately of Drude form 



(20) 



^ Let us consider limr^s ri-Vf G(r, s')a(s') ds' for a coordinate 
system where ri = e^, r = (0, 0, 2;), and s' = p(cos 0, sin 0) 
is given in polar coordinates p and (j). We compute the bound- 
ary integral within a small circle with radius within which 
the surface charge cr can be approximated by a constant. The 
integral then becomes 

/r — s' 3 
— ds' lim 27TZ I pdp (p^ + z^)" 2 = 27t . 
\r- s'\^ z^o Jo 



^ For the diagonal elements of G we change to polar coordinates p 
and (f), and perform the integration within a triangle where r((f)) 
denotes the upper limit of p for a given cf). We then obtain 

/ d(f) pdp — = dcf) , 

Jo Jo p Jo 

which can be easily calculated numerically. The diagonal ele- 
ments of F are all zero. 
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Here eo is a dielectric constant, accounting for screening 
of d-band electrons in transition metals, cOp is the bulk 
plasma frequency, and 7 the Landau damping of plas- 
mons. In gold the above Drude form is valid for energies 
below the threshold value of approximately 2 eV where 
d-band excitations set in. Typical values for Au, which 
we will be using in this work, are eo = 10, cjp = 9 eV, 
and 1/7 = 10 fs. Instead of the Drude form we can also 
use a hydrodynamic model where electrons with particle 
density no move freely in a medium with the background 
dielectric constant eo- This model is completely equiv- 
alent to the framework of the Drude dielectric function. 
As we will show below, through the hydrodynamic model 
we can establish a microscopic description for the elec- 
tron dynamics, which will prove useful for interpreting 
our results in physical terms and for performing a quan- 
tization of surface plasmons. 

The energy of a classical electron plasma is the sum of 
kinetic and electrostatic energy [27l [28] 



(21) 



Here p{r) is the charge density displacement from equilib- 
rium, ^{r) is the electrostatic potential induced by 
and ^(r) is the velocity potential, whose derivative gives 
the velocity density v = — [28]. For the surface plas- 
mons of our present concern we consider surface charge 
distributions a which are nonzero only at the surface of 
the MNP. As detailed in appendix [Dj the Hamilton func- 
tion (21 ) can be rewritten in a boundary element method 



approach as 



H 



.2 



27r(eo + e6)l + (eo-e6)F 



Go- 



}. (22) 



Here a is the vector of the surface charges within the 
discretized surface elements, cOp = {Atttio)^ is the plasma 
frequency, and F is the surface derivative of the Green 
function with respect to the second argument. Equa- 
tion (22) accounts for undamped plasma oscillations. 



Damping as well as coupling to the molecular dipole has 
to be included additionally, as we will discuss further be- 
low. 



Before doing so, we show that Eq. (22) allows us to 



compute eigenmodes for the surface charge oscillations. 
As the matrices appearing in the Hamilton function are 
real and symmetric, they can be diagonalized simultane- 
ously. Let cj^ and ux denote the eigenvectors and eigen- 
values of the generalized eigenvalue problem 



(^27rl + 



27r(eo + e5)l + (eo-e5)F Gux- (23) 



The eigenvectors ux can be chosen real and are orthog- 
onal in the sense vy^ {2it1^F)-^Gux' = /3a ^^aa', where 



Px is a real number. We can now expand the surface 
charge distribution in terms of these eigenfunctions viz. 
[271 [281 [38] 



/ 2no 
V^a/^a 



^a(^) {; 



-iojxt 



(24) 



where ax are expansion coefficients for the plasmon mode 
A. With Eq. (24) the Hamilton function can be brought 
to the intriguing form H = ^ uJx{cL\a^ + ^x^x)- "^^^ 
quantization may now be performed through the substi- 



tution a^^ a\ 



where a^ and at are the oper 



ators for annihilation and creation of a surface plasmon 
mode A, respectively. They obey the usual bosonic com- 
mutation relations. With these quantized surface plas- 
mon modes Eq. ([23]) can be brought to the form 



(25) 



Hpi = ^u;x (4^A + i) • 



This is the Hamilton operator for surface plasmons. In 
principle it can be used to describe the non-linear prop- 
erties of metallic nanoparticles, although in this work we 
shall only be interested in linear response. What still re- 
mains to be done is to show how Landau damping can 
be included, and how to obtain the nonradiative and ra- 



diative scattering rates starting from Eq. (25). 

Let us first consider the coupling between the molec- 
ular dipole and the surface plasmons (24). Using the 



expressions given in appendix |D] we can compute for a 
given surface mode ux the scalar potential $a and its 
surface derivative <I>^ according to 



^1 = 



-4n 
Q- 



27r(eo + eb)l + (eo-eb)F 
1 {2iTt - f) $A • 



(26) 



The scalar potential <I>a at the position of the molecu- 
lar dipole is then obtained by means of Green's second 
identity [Eq. (Dl)], and the coupling between the molec- 
ular dipole and the surface plasmon ux from the standard 
expression gx = Pdip^^A- In our numerical approach we 
again describe the charge distribution pdip of the molec- 
ular dipole by means of two oppositely charged particles 
separated by a small distance. Within this approach we 
then obtain for the coupling between the molecule and 
the surface plasmons the interaction Hamiltonian 



pi — mol 



A 



(^A^+^A +^A^-4) ' (27) 



where we have made use of the rotating- wave approxima- 
tion. d-± are the molecular excitation and de-excitation 
operators introduced in Eq. (T]). The first term on the 
right-hand side of Eq. ( [27| ) describes the creation of a 
surface plasmon through de-excitation of the molecule, 
and the second term the reversed process of plasmon an- 
nihilation and molecule excitation. 
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FIG. 1: (a) Comparison of BEM and Mie calculations for 
scattered light intensity of a spherical gold nanoparticle, as 
computed within the quasistatic approximation and for the 
dielectric function of Ref . [39 . The symbols represent results 
for different surface discretizations with (c) 144 vertices, (d) 
400 vertices, (e) 400 vertices and a finer grid around the north 
pole, and (f) 900 vertices, (b) Nonradiative decay rate of a 
molecule in vicinity of the spherical nanoparticle, in units of 
free-space decay rate 7° and for a molecular dipole oriented 
in radial direction. The transition energy of the molecule 
is 2.4 eV [see arrow in panel (a)], and the diameter of the 
nanosphere is 10 nm. The inset reports the ratio of the BEM 
scattering rate with the exact Mie result. 



2. Mie theory 

For spherical particles it becomes possible to obtain an- 
alytic expressions for the surface plasmon energies (jJ\ and 
coupling constants g\ by means of Mie theory. We here 
provide for comparison the corresponding results which 
were first obtained by Ritchie [27 . The plasmon energies 
read 



{I + 1)65 + leo ' 



(28) 



where / denotes the usual degrees of spherical harmonics. 
For the coupling constants we get 




n}yo}^^^ air, 



(29) 



Here a is the radius of the sphere, and aim are the expan- 
sion coefficients of the molecular dipole given in appendix 

o 



3. Surface plasmon dynamics 

With the quantization of the surface plasmons we have 
opened the quantum optics toolbox, which would allow 
us to investigate strong coupling or other nonlinear ef- 
fects |33]. In this work, however, we will use the surface 
plasmon eigenmodes only to analyze the behavior of the 
radiative and non-radiative decay rates 7^ and 7nr- Pur- 
suing a similar open-system approach as in appendix [A} 
but for the molecule-plasmon coupling (27), we can ex- 
press the nonradiative decay rate as 



E 



\9\\ 



-iujt 



a^{{))a\{t)) dt. (30) 



Here the exponential accounts for the free propagation of 
the molecule, and the terms in brackets describe the plas- 
mon correlation function to be evaluated at zero temper- 
ature. The latter is of the form exp[i(a;A — ^2)^]' where 
7 is the Landau damping of the plasmons. In a more 
general approach we could have also introduced damp- 
ing through a master equation approach with suitable 
Lindblad operators. We can finally solve the scattering 
integral of Eq. (ISOl) to arrive at 
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(31) 



Again we have neglected the frequency renormalization 
due to plasmon coupling. Eq. (31) shows that, through 



an eigenmode expansion, the nonradiative decay rate can 
be decomposed into the different contributions for each 
plasmon mode. Each component is given by the square of 
the coupling constant together with a Lorentzian whose 
broadening is given by the Landau damping. Note that 
in the limit 7^0 the Lorentzian would reduce to 
t^5{(jO\ — uo). In addition to Eq. (31), the radiative scat- 



tering rate can be obtained by computing, within lowest- 
order perturbation theory, the occupation of the plasmon 
modes, and multiplying it with the dipole moment of the 
corresponding mode. 



III. RESULTS 
A. Spherical particles 

Figure [TJa) reports the scattered light intensity from 
a spherical gold nanoparticle, as computed within our 
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FIG. 2: Scattered light intensity for different nanoparticle 
shapes of (a) rod, (b) cigar, (c) elUpsoid, (d) disk, and (e) 
sphere. The ratio of height : diameter is 5 : 1 for (a,b) and 
1 : 5 for (c,d). The soUd and dashed hnes show results ob- 



tained for the dielectric function 39 and the Drude form (20 ), 



respectively. Light polarization is chosen along the long axes. 
The dots above the spectra indicate the positions of the sur- 
face plasmon eigenmodes, as obtained from the solutions of 
Eq. ([23|. The corresponding eigenmodes of lowest energy are 
shown at the bottom of the figure. For the considered par- 
ticle shapes only the eigenmodes with lowest energy have a 
non-zero dipole moment and can couple to light. 



quasistatic approach, and for the dielectric function of 
Ref. [39^. The solid line is the result of Mie theory, 
which is summarized in appendix [Cj and the symbols 
show results of our boundary element method approach 
described in Sec. |IIC[ The broad peak at a photon energy 
of approximately 2.3 eV is associated to the dipole res- 
onance of the nanoparticle. We observe from the figure 
that the difference between BEM and Mie results reduces 
with increasing number of surface elements. 

Panel (b) shows the nonradiative decay rate for a 
nanosphere, with 10 nm diameter, and a molecule placed 
at a given distance away from the particle. We assume 
that the dipole is oriented in radial direction, and set 
the molecular transition frequency to 2.4 eV. One ob- 
serves that with decreasing distance the decay enhance- 
ment drastically increases, reaching a value of approxi- 
mately 10^ for the smallest distance of 0.2 nm. The sym- 
bols are the results of our BEM calculations. From the 
inset, which reports the ratio between the BEM and Mie 
decay rates, one observes that the BEM results give reli- 
able results only down to molecule-nanoparticle distances 
comparable to the discretization length. The results for a 



sphere discretization with 144 vertices [see panel (c)] sig- 
nificantly deviates already at a distance of 2 nm, while 
the results for sphere discretizations with 400 and 900 
vertices [panels (d) and (e)] show nice agreement down 
to values of about 1 nm. Things improve considerably 
for the special sphere discretization shown in panel (e). 
Here we have introduced a finer mesh around the north 
pole where the molecule approaches the sphere. As evi- 
dent from the inset of panel (b), the difference between 
the BEM and Mie calculations is at most ten percent, 
even for such small distances as 0.2 nm. 



B. Other particle shapes 

In the following we consider different particle shapes 
where the scattering properties cannot be obtained an- 
alytically. Figure [2] reports optical spectra for different 
particle shapes, which are shown in the lower panel. In all 
cases the light polarization is along the long axes of the 
nanoparticles. The plasmon energy of the dipole mode in- 
creases from the quasi one-dimensional rod (a) and cigar 
(b), over the quasi two-dimensional ellipsoid (c) and disk 
(d), to the sphere (e). The solid and dashed lines show re- 
sults as obtained from the dielectric function of Ref. [39] 
and the Drude form ( po] ), respectively. With the excep- 
tion of the sphere, both results are in nice agreement and 
thus justify the Drude description for quasi one- and two- 
dimensional nanoparticles. From the comparison of the 
results for the rod and cigar, as well as for the ellipsoid 
and disk, we observe that the detailed shape of the par- 
ticle has no dramatic impact on the spectra. Finally, the 
height of the peaks has the approximate ratio 1 : 4 : 9 for 
the particles of different dimension. Indeed, this is the 
behavior one would expect for a dipole radiator where 
the oscillator strength is distributed between (a,b) one, 
(c,d) two, and (e) three dipole modes, where only one is 
optically excited. 

The dots in the upper panel indicate the plasmon 
eigenenergies ujx obtained from the generalized eigen- 
value problem (23). The lowest dipole modes are pre- 



cisely at the positions of the corresponding peaks in the 
optical spectra. In the lower panel of the figure we show 
the eigenmodes ux for the two surface plasmons of lowest 
energy. For the particle shapes under consideration only 
the mode of lowest energy has a finite dipole moment 
and can thus couple to light. A striking feature of the 
eigenenergies is that with decreasing energy of the dipole 
mode the energy separation to the next mode increases. 
It is of the order of 0.4 eV for the quasi one dimensional 
particles (a,b), about 0.2 eV for the quasi two dimen- 
sional particles (c,d), and decreases by a further factor of 
two for the sphere. For the sphere we also observe a quasi 
continuum of surface plasmon modes at energies around 
2.5 eV, in agreement with the result (28) of Mie theory. 



As we will show next, these different energy separations 
between plasmon states have important consequences for 
the non-radiative decay rates. 
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FIG. 3: (a) Non-radiative decay rate in units of radiative 
free-space decay rate 7° for different nanoparticle shapes of 
sphere (circles), disk (triangles), and cigar (squares). The 
length of the particles along the long axes are 28 nm for the 
cigar, 16 nm for the disk, and 10 nm for the sphere. These 
dimensions have been chosen such that all particles have the 
same volume. The orientation of the molecular dipole and the 
direction along which the dipole approaches the nanoparticle 
are indicated in the inset, and the molecular transition energy 
is assumed to be in resonance with the dipole modes, (b) 
Quantum yield ^r/lnr for the particles shown in panel (a), 
and for different molecule-nanoparticle distances measured in 
units of the sphere diameter of 10 nm. (c) Decomposition of 
the non-radiative decay rate into eigenmodes, as described in 
text. 




FIG. 4: Enhancement of FRET transfer rate for two differ- 
ent donor positions. The transition frequency of the donor 
molecule is in resonance with the dipole surface plasmon. We 
perform an average over all dipole orientations of the donor 
and acceptor molecules. Note that in the quasistatic limit the 
FRET enhancement does not depend on the particle size. 



tor of, e.g. 2 or 4, ^nr would decrease by a factor of 8 and 
64, respectively, and the quantum yield increase by the 
same factor. Thus, for a quasi one- or two-dimensional 
metallic nanoparticle whose length is around 50 nm, the 
quantum yield can easily be of the order of several ten 
percent. 



We can now use Eq. (31) to decompose the scattering 
rate into contributions for different plasmon modes A. 
This is what is shown in the lower panels of the figure for 
(c) the cigar, (d) the disk, and (e) the sphere. Each point 
in the density plots corresponds to a given molecule- 
particle distance and plasmon energy, and we have in- 
troduced a broadening of the order of 7 for each plasmon 
mode. The colors indicate the relative importance of the 
different plasmon modes. For the cigar-shaped particle, 
we observe in panel (c) that up to very small distances it 
is primarily the dipole surface plasmon mode that couples 
to the molecule. This is because of the large energy sepa- 
ration between the plasmon modes of lowest energy, and 
the resulting weak coupling to the off-resonant modes. 
Only for the smallest distances the coupling constants g\ 
dominate over the detunings uox — uo. Similar behavior is 
observed in panel (d) for the disk-shaped particle. Things 
are different in panel (e) showing results for the sphere. 
Because of the small energy separation between the dif- 
ferent surface plasmon modes, already at relatively small 
distances, say around 0.1 in units of the sphere diameter, 
the molecule starts to significantly excite the higher-lying 
plasmon modes. As consequence, the quantum yield dra- 
matically drops at these distances. 



C. Forster resonant enery transfer (FRET) 



Figure [3] shows (a) the non-radiative decay rate, and 
(b) the quantum yield "^r/lnn foi" the different particle 
shapes shown in panel (a). In all cases the molecular 
transition frequency is assumed to be in resonance with 
the dipole modes. Note that the volume of the particles 
corresponds to that of a relatively small sphere of 10 nm 
diameter. By increasing the size of the particles by a fac- 



We conclude this section by briefly discussing Forster 
energy resonance transfer (FRET) in presence of 
nanoparticles. Here a donor and acceptor molecule ex- 
change energy through dipole-dipole interaction. This 
process become drastically enhanced if the molecules can 
benefit from the nearfield enhancement of the metal- 
lic nanoparticle jlQl |41|. For given donor and acceptor 
dipole moments /L^don and /L^acc, the enhancement of the 
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FRET process can be computed from 



APPENDIX A: MASTER EQUATION 



^FRET 



I Mace ' ^{'^acc-) ''^don) ' Mdoi 
iMaee ' ^o(''^aee7 ''^don) ' Mdo 



|2 • 



(32) 



rdon and race denote the positions of the donor and ac- 
ceptor molecule, and G and Go are the dyadic Green 
functions in presence and absence of the nanoparticle. In 
our approach, they are computed in the quasistatic ap- 
proximation and for the molecular transition frequency, 
in accordance to the prescription of Sec. II C[ Note that 
the denominator in Eq. (32) is precisely the free-space 
dipole-dipole coupling. 

Figure |4] shows the FRET enhancement for a disk- 
shaped nanoparticle and for a molecular transition fre- 
quency equal to the dipole surface plasmon energy. The 
positions of the donor molecules are indicated in panels 
(a) and (b). For simplicity, in the plot we assume that 
the acceptor molecules are placed on a sheet, and we have 
averaged over all orientation of the donor and acceptor 
molecules. The FRET transfer becomes enhanced by up 
to four orders of magnitude, in particular for transfer 
processes across the nanoparticle. Here, the donor and 
acceptor molecule exchange their energy via the resonant 
surface plasmon mode. Note that in order to estimate the 
FRET efficiency we additionally have to consider the en- 
hancement of the radiative and nonradiative decay pro- 
cesses, as discussed in Ref. j4lj . 



IV. SUMMARY 

In summary, we have developed a general framework, 
based on a suitable description of quantum electrody- 
namics in the presence of absorbing media and a general 
boundary element method approach, for the description 
of the modified molecular dynamics in presence of metal- 
lic nanoparticles. Within the quasistatic limit the decay 
rate has been separated into radiative and non-radiative 
decay channels, which both obey a simple scaling behav- 
ior. We have discussed different nanoparticle discretiza- 
tions, and have shown that reliable results can only be ob- 
tained for properly defined surface discretizations. From 
a comparison with Mie theory we have estimated an er- 
rorbar of less than ten percent for our numerical results. 
We have performed an eigenmode expansion for the sur- 
face plasmons, and have shown that the higher quantum 
yield for quasi one- and two-dimensional nanoparticles 
can be attributed to the large energy splitting of the cor- 
responding eigenmodes. Finally, we have shown that sig- 
nificant enhancements of the Forster transfer rates can 
be expected for donor and acceptor molecules placed in 
the vicinity of the nanoparticle. 



Acknowledgments 



In this appendix we provide some details for the deriva- 
tion of the Fermi-golden-rule result ([3|. As we are dealing 
with an open system, i.e., a molecule interacting with the 
continuum of photon modes, we have to adopt a density- 
matrix description [33^. We assume that the molecule is 
initially in the state 1. The corresponding density opera- 
tor is of the form po = |1)(1|. The total density operator 
w has to additionally account for the degrees of freedom 
of the electromagnetic environment. Let us make the 
usual assumption that at time zero the molecule is de- 
coupled from the electromagnetic environment. The total 
density operator then factorizes wq = po ^ pji^ where pR 
is the density operator of the electromagnetic reservoir 
which we consider to be in thermal equilibrium. The time 
evolution of w{t) follows from the Liouville von-Neumann 
equation w{t) = —i[Hop{t)^w{t)]^ which we have writ- 
ten in the interaction representation with Hop the light- 
matter coupling of Eq. ([2|. As we are only interested 
in the molecular dynamics, we remove the explicit de- 
pendence on the reservoir by tracing out the degrees of 
freedom of the electromagnetic environment through 

^{t) = tTRiij{t) = -itiR [[Hop{t),w{t)]^ . (Al) 

So far we have not achieved too much since the expression 
on the right-hand side still involves the full density opera- 
tor w. To proceed, we describe the light-matter coupling 
in lowest order time-dependent perturbation theory, also 
known as the Born-Markov approximation [33^, 



m 



Jo 



Hop{t), Hop(t'),p(t)® PR 



dt'. 



(A2) 

Equation (lA2| is correct up to second order in the light- 
matter coupling Hop. There is a subtle point regard- 
ing the time argument of p{t) on the right-hand side, 
for which we refer the interested reader to the literature 
[42, 43 . We next insert the light-matt er co upling ([2| in 
rotating- wave approximation into Eq. (A2), to arrive at 
the generalized master equation 



m 



dt' 

Et{t)E-{t'))d-{t)d+{t')m 

Et{t')E-{t)) md-{t')d+{t) 

E+{t')E-{t))d+{t)md-{t') 
Et{t)E-{f))d+{t')p{t)d-{t) 

(±) ^ (T)l . (A3) 



We are grateful to Joachim Krenn and Alfred Leitner 
for most helpful discussions. 



Here we have introduced the shorthand notation ( . ) = 
tri? {pR . ), have exploited the cyclic permutation of the 
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field operators under the trace, and have used that the ex- 



pectation values {E- E- ) vanish in thermal equilibrium. 
The four terms explicitly given on the right-hand side 
of Eq. (A3) correspond to photon emissions, and those 



resulting from the exchange of (±) with (t) to photon 
absorptions. As in this work we will be dealing with elec- 
tromagnetic fields in the visible regime and room temper- 
ature, we can safely ignore the thermal occupation of the 
radiation modes, and thus neglect absorption processes 
of thermal photons throughout. 

Let us assume that the dipole operators oscillate 
with the molecular transition frequency uoq according to 
d^{t) = e^*^o^d±. In addition, we introduce the Fourier 
transform of the field operators 



poo 

Jo 



2^ 



(A4) 



We will use the same symbol E for the field operators in 
time and frequency space as long as there is no danger of 
confusion. For sufficiently large times t, where uo^t ^ 1 
is fulfilled, we can replace the lower limit in the time 
integration (A3) by — oo. This replacement is called the 
adiabatic approximation [33l|43] and is usually considered 
to be valid under broad conditions. The time integrations 
then reduce to integrals over exponentials, which can be 
solved, with the relative time t = t — t' ^ according to 

^0 

g-i(a;o-a;+iO+)r^^ = ^ ^ 1t5(uJ - CJq) • 

CJO - ^ + ^0+ ^ ^ 

(A5) 

The infinitesimally small and positive quantity 0+ has 
been introduced to provide a damping of the exponen- 
tial at early times. Its appearance has nothing to do 
with a real damping process. Rather it allows to per- 
form the adiabatic limit in the integration prior to the 
thermodynamic limit in the expectation value (.) [44] . 
Finally, in the last expression on the right hand side of 



Eq. (A5) we have only considered Dirac's delta function, 
which accounts for the energy conservation in the scat- 
tering processes, and have neglected the Cauchy princi- 
pal value. The latter contribution is responsible for an 
energy shift of the molecular transition, usually known 
as the Lamb shift [33], which could be incorporated into 
a renormalized ujq. In thermal equilibrium, the expec- 
tation value of the field operators {Ef {ijo)EJ {uo')) = 

2ii8{uj—u'){E^{uj)Ej{uj)) is diagonal in frequency space. 
Putting together all the results and changing back from 
the interaction to the Schrodinger representation, then 
brings us to the master equation 



^{Et{uo)E-{u^) 



d. dj p^pd~^ 



2d+ pdi 



(A6) 



The first term on the right hand side accounts for the 
free propagation of the molecule, and the second term 
for scattering processes due to the coupling to the elec- 
tromagnetic fields. The terms in parentheses are remi- 
niscent of the Lindblad master equation. Finally, using 



the explicit form ([T]) for the dipole operators and tak- 
ing the matrix elements pn = of the operator 



equation (A6), we find that the upper-state population 
of the moL 



•lecule decays according to pn 
the decay rate 7 given in Eq. (|3|. 



-7Pii 



with 



APPENDIX B: CURRENT CORRELATION 
FUNCTION 

In this appendix we provide some details of how to re- 
late the current correlation function {j^ {r^uj)jf {r' ^uj)) 
to the imaginary part of the dielectric function of the 
metal. We will only consider local and isotropic me- 
dia, such that the expression given above reduces to 
^ki^{T^ — '^'){j^{<^)j~ {'^))' The calculation of such cor- 
relation functions is a common problem in solid state 
physics j45ll46] . 

We start with a general result of linear response theory 



/// N 47r 
e iuo) = -^\sm 



f 

Jo 



e'^U[j+{t),j-m)dt, (Bl) 



which relates the imaginary part of the dielectric func- 
tion to the retarded current correlation function j34] . 
A common link between the ordered and retarded cor- 
relation function is provided by the spectral function 

p{t) = ([j {t) ^ (0)]^ . To compute p{t) we would need a 

microscopic model for the electron dynamics in the metal, 
which is certainly beyond the scope of our work. How- 
ever, in order to obtain, in thermal equilibrium, the rela- 
tion between the different correlation functions it suffices 
to assume that the states |m) and energies Em of the 
full many-body hamiltonian H could be determined in 
principle. In fact, we never have to compute these states 
explicitly. Inserting in ( . ) = tr(e~^^ . ) a complete 
set of states, where /3 = l/{kBT) is the inverse tempera- 
ture and Z the partition function, we then obtain for the 
Fourier transform of the spectral function (45] [46] 



(B2) 



pM = (l-e-''-)Z-i^e-''^'"|HMn) 

m,n 

x2TT5{u-[En-Ej^). 



Prom Eq. (B2) one then obtains {j (w)) 



p{uj)/{l — e~^"). Similarly, upon insertion of the many- 
body states into the retarded correlation function we get 
p{uj) = u!^ e" (uj) / {2Tr) . Putting together all results we 
arrive at our final expression 



27r ' 



(B3) 



with nth(^) = l/(e^^ — 1) the Bose-Einstein distribution 
function. For the problem of our present concern we can 
safely neglect the thermal occupation nth of photons. 
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APPENDIX C: MIE THEORY IN THE 
QUASISTATIC LIMIT 



plasma. Let us consider first the potential energy term. 
Our starting point is provided by Green's second theorem 



We consider the situation where a single molecule with 
dipole moment /j. is located at some distance away from 
a spherical nanoparticle with radius a. We expand the 
potential inside and outside the particle in terms of spher- 
ical harmonics Yi^ [35^ 



OO I 

^=0 m=-l 
CO I 



47r 



2. 



47r 



^ 21 

1^0 m^-l 



1 



Clr 



^dip . (CI) 



Here bim and cim are expansion coefficients, to be de- 
termined below, r is the distance from the center of the 
sphere, and <l>dip is the potential of the molecular dipole. 
^dip can be expanded similar to the first expression for 
^out5 but with different expansion coefficients aim- They 
read [47] 



+zv'KlTl)(nXM)-Xr^(^,0)}. (C2) 

Here n is the unit vector from the center of the sphere to 
the dipole located at the position characterized by r, 6>, 
^, and Xim{0^ (j)) are the vector spherical harmonics [35 . 
We can now use the boundary conditions of Maxwell's 
theory to relate the different coefficients viz. 



21+1 



(1 — er)la 
(l + e^)/ + l 



(C3) 



with €r = €m/^b being the ratio between the dielec- 
tric functions of the MNP and the background material 
within which the particle is embedded. The electric field, 
needed for the calculation of the non-radiative decay rate 
( 11 ), becomes 



^=0 m=-l 



47r Clr, 



21 + 1 



{{l^l)flYirn{e,^) 



+ 1) n X Xl^{e, (/))} + E^^ . (C4) 

For the dipole moment of the MNP, induced by the 
molecule, we get [35 



Mmnp 



Here is the unit vector along z, and e± 
iey)/^/2. 



(C5) 



47rcE>(r) = / {G{r,s') 
Jon ^ 



d^{s') dG{r,s') 



dh' 



dh' 



^{s')]ds' 



(Dl) 

where the partial derivatives denote the surface deriva- 
tives h • V w ith respect to s' . Simil arly to the discussion 
in Sec. II we can perform in Eq. (Dl) the limit r ^ s 
to obtain a relation between <1> and its surface deriva- 
tive. Within our boundary element method approach we 
obtain (27rl ± F)<1> = ±G^' for the hmit taken in- or 
outside the metallic nanoparticle. F denotes the surface 
derivative of G with respect to the second argument. To 
obtain a relation between $ and the surface charge a, we 
make use of the boundary condition n-{Di) — Dm) = 47r<j. 
Together with Eq. (Dl) we then obtain 



$ = 4^{2^(6o + e5)l + (eo-65)F} ^Ga. (D2) 

As for the kinetic energy of the classical plasma, we use 
the continuity equation 



dtTi — — noV • V — noV^^ , 



(D3) 



to relate the density displacement n and the velocity po- 
tential ^. Integration of the continuity equation (D3) 
over a small cylinder Q (height h ^ and base 5S) con- 
taining a small surface element, then gives for the right 
hand side of Eq. ([D3l 



/ noV^^dV= / non-V^ds^no^6S. (D4) 
Jn Jon dn 



Here, |¥ = n • denotes the surface derivative of the 
velocity potential. Together with the left-hand side of 
Eq. (D3) we find the link between ^ and a. 



no 



dh 



(D5) 



Using Green's second identity, we can again relate within 
our boundary element method approach the surface 
derivative of ^ to the velocity potential viz. ^ = 
(27rl+i^)~^G ^^ We can thus express the kinetic energy 



L 



F)-^GcT. (D6) 



APPENDIX D: HYDRODYNAMIC MODEL 



In this appendix we show how to obtain Eq. (22) start- 



through the time derivative of the surface charge distri- 
bution. Putting finally together the contributions for the 
kinetic and potential energies, we arrive at expression 
( 22 ) for the energy of a classical plasma in terms of the 



ing from expression (21) for the energy of a classical surface charge a. 
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